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ABSTRACT 


Rough surface scattering plays a crucial role in the statistics of acoustic propagation 
signals, especially at mid-frequencies and higher (e.g., acoustic communications 
systems). For many years, the effects of rough surface scattering were computed using 
simple models that were applied in two dimensions (2-D) only. 

A prescribed method of computing 2-D rough surface scattering directly in a 
parabolic equation model based on the Split-Step Fourier algorithm was introduced by 
Tappert and Nghiem-Phu in the mid-1980s. This method has been successfully 
implemented in various 2-D parabolic equation models, including the Monterey Miami 
Parabolic Equation model. However, some scientific research of more formal scattering 
predictions have suggested that out-of-plane, three dimensional (3-D) scattering may lead 
to significant disparities in the scattered field statistics. Introducing a hybrid 
implementation for the scattering effect in the field transformation equations using a tri¬ 
diagonal solution with the Pade approximant to obtain a system of equations for 
azimuthal corrections will support predictions of the effect of surface scattering on 3-D 
propagation, which is critical in evaluating the variability in underwater acoustic 
propagation. Results of the 3-D scattering calculations obtained are compared with the 
output of basic 2-D interface perturbations utilizing the standard 2-D approach. 
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I. 


INTRODUCTION 


In reality, all surfaces are rough, and it is often important to know how surface 
roughness changes the properties of scattered acoustic waves (SAW) compared with 
those reflected from a plane surface. The extent to which surface roughness affects wave 
scattering and finally the way it propagates is therefore of great interest. Over the last few 
decades this problem has received considerable attention in several branches of physics, 
meteorology and engineering. Of major interest has been the study of electromagnetic 
wave scattering from rough surfaces, which has been crucial in designing many different 
operational systems (e.g., radar systems) where the effects of the roughness of the land 
and sea surfaces on signals have been studied (e.g., Ogilvy 1987). Similarly, considerable 
interest has been concentrated on the interaction of sound waves with rough surfaces 
owing to its significance in several fields. Rough-surface effects on the acoustics field 
have been extensively studied as well (e.g., Chimenti and Lobmis 1998). 

Acoustics, generally speaking, is one of the primary means by which the 
underwater environment can be sampled, because acoustic waves are used in the 
underwater medium in much the same way that electromagnetic waves are used in the 
atmosphere. The propagation of acoustic waves in the underwater environment is 
important for several areas, including navigation, ocean observation, detection of 
submerged objects, and real-time underwater communications. As the ocean environment 
changes, the refraction and scattering of the acoustic field is affected. Over long time and 
space scales, the change in the characteristics of the sea-water throughout the water 
column causes the sound speed profile to change. This in turn shifts the acoustic paths 
and changes the overall capacity for the ocean to transmit acoustic waves. 

Modulating a pressure-release ocean surface over short time and range scales has 
been found to be more efficient in studying the acoustic paths in which virtually all 
incident acoustic energy is reflected in the water column (Medwin and Clay 1998). 
Consequently, multi-paths will be generated as long as the surface modulates over time 
and space (Senne 2012). 
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Certainly, there are numerous environmental factors that can influence the 
acoustic wave propagation in the ocean, such as water temperature and salinity, as well as 
sea bottom and surface perturbations. The impact of all these factors on acoustic wave 
propagation are studied in order to determine the performance level of underwater 
acoustic systems and also to develop techniques for using acoustic signals to measure the 
physical parameters of the ocean. Surface roughness is just one among several 
environmental parameters that can have noticeable and essential influence on underwater 
acoustic propagation. 

Biot’s 1957 paper suggested one way to approach the study of the scattering of 
acoustic waves from a rough surface. The idea was to treat the surface as a collection of 
hemispherical protuberances—or “bosses.” He states: 

A simple solution is developed for the reflected waves on a rough surface 
from a simple hannonic point source. It is assumed that the roughness is 
represented by a distribution of hemispherical bosses whose size and 
mutual distance are small relative to the wavelength. It is shown that under 
these conditions the effect of the roughness is equivalent to a boundary 
condition for the wave equation. This boundary condition embodies the 
surface polarization and the mutual interaction of the bosses. 

Until recently, only 2-D models have been implemented in order to compute and 
illustrate the rough surface effects upon the acoustic radials within the field (e.g., 
Shchegrov 1995). More recent 3-D modeling of rough surface scattering has been based 
on ray models or localized calculations of full-field scattering. As full field 
approximations provide more information than ray path solutions, a parabolic equation 
(PE) approximation is chosen to model the acoustic field out to long range with many 
multipath contributions from bottom and surface scattering. 

In this study, a three dimensional (3-D) ocean acoustic model has been combined 
with a 2-D rough surface fonnulation of a parabolic equation model to simulate the 
effects of sea surface roughness on acoustic wave propagation with multiple surface 
interactions. Rough sea surface realizations are generated in the acoustic model. An 
approach to achieve high resolution and accurate results while maintaining computational 
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efficiency is applied. The results of the scattering calculations are then compared with 
results from the simple 2-D model. 

Throughout this thesis, the PE model is used to detennine the acoustic field at a 
variety of range and depth locations in an underwater environment. The second chapter 
reviews the past research and the analytical tools developed. It also includes the 
derivation of the PE equations for the 3-D scattering model approach and a description of 
its usage - including the treatment of rough sea surfaces. Concerning the numerical 
implementations, Chapter III introduces the hybrid implementation for the scattering term 
using a tri-diagonal solution in addition to defining the rough surface perturbation 
equations that are required for the solution. Results from this acoustic model are 
presented in Chapter IV, which includes an accounting of the field study through 
different surface perturbations by introducing five separate models while analyzing the 
effects of rough surface on acoustic propagation for each model. The PE model coupled 
with each type of 3-D scattering model is compared against previous 2-D model results in 
Chapter IV. 

A. PROBLEM STATEMENT 

The intensity range of the reflected waves can vary widely due to the impact of 
the evolving sea surface over a short time scale (Ogilvy 1987). Generally, acoustic waves 
scatter from surfaces with wavelengths that are comparable to or larger than the acoustic 
wavelength. As the acoustic wave frequency increases, the associated acoustic 
wavelength decreases. Consequently, smaller scales of roughness at the boundaries 
become more important. This is a key point when considering sea surface modeling. Sea 
surfaces are made up of a variety of wave amplitudes and wavelengths that are 
superimposed. Depending on the acoustic frequency being used, several of these 
individual surface wave components may play important roles in affecting the acoustic 
wave. For example, high frequencies, upon which acoustic communication systems 
depend, are more susceptible to small fluctuations in the sea surface roughness (Senne 
2012 ). 
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The importance of the 3-D out-of-plane models has been suggested by Karjadi 
(2011), as he states: 

Although the present 2-D acoustic model provides a simple and 
computationally efficient technique of examining the fluctuations of 
arrival times and angles caused by surface waves, the 2-D approach does 
not directly address the out-of-plane scattering caused by surface waves; 
to properly model these out-of-plane scattering effects, a full 3-D acoustic 
wave model is required. 

Moreover, Senne (2012) has discussed the importance of the realistic evolving 
rough sea surface by stating: 

Understanding acoustic interactions in the vicinity of the rough sea surface 
and maximizing the usability of underwater communication systems can 
be achieved through the modeling of acoustic propagation within realistic 
ocean environments that include time-varying rough sea surfaces. 

B. PURPOSE STATEMENT 

The goal of the present work is to realize the general effects of the environment 
variability, such as rough surface scatter, on the 3-D propagation which is critical in 
evaluating the limitations of existing sonar prediction models. If the 3-D effects are found 
to be significant, this suggests sonar prediction models must be adapted to account for 
these effects. 

In Chapter II, the selected 2-D and 3-D techniques for rough surface scattering 
calculations are introduced along with a description of how the two models are 
integrated. 

Comparisons between this model output and the output of basic 2-D interface 
perturbation models are considered to show the feasibility of the approach and 
demonstrate the ability of the model to track acoustic wave propagation under fluctuating 
sea surfaces from one range step to another. 
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II. THEORY 


A theoretical procedure using a full 3-D PE model is developed in an effort to 
solve the acoustic propagation problem at the surface of the ocean in the presence of 
surface displacements. This theory is believed to be an improvement over past efforts that 
computed the effects of rough surface scattering using simple models that were applied in 
2-D only. 

A. REVIEW OF PAST RESEARCH AND ANALYTICAL TOOLS 

DEVELOPED 

We discuss briefly here each of the major techniques which have been applied to 
the issue of rough surface scattering. Due to the subject importance to the various 
scientific fields, there are numerous papers related to rough surface scattering. We have 
selected a representative sample in order to introduce previous techniques used and the 
analytical tools developed. 

Probably the first reported work on scattering from rough surfaces was by 
Rayleigh (1878) who discussed the scattering of a normally incident sound wave by a 
corrugated (sinusoidal) surface separating two acoustic media. He assumed that the 
scattered field could be written as a sum of plane waves travelling away from the rough 
surface. This approach has since been extended by other authors to non-normal incidence 
and to random rough surfaces (Ogilvy 1987). 

In the perturbation approach, the rough surface is regarded as a perturbation to a 
smooth plane, and the consequent change in the scattering coefficient which measures 
how much of the incident acoustic waves are scattered, due to the presence of roughness, 
is calculated. This approach requires that the height deviation of the surface away from 
the smooth plane is everywhere small compared to the wavelength of the incident wave. 
In addition, the slope of the surface must be small in comparison to unity (Gilbert and 
Knopoff 1960). 
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Undoubtedly the approach most prevalent in the literature is the tangent plane, or 
Kirchhoff, method. Probably the first use of this theory was by Eckart (1953), applied to 
acoustic wave scattering, and Davies (1955) for electromagnetic wave scattering. Here 
the scattering surface is assumed to be everywhere sufficiently smooth so that the 
reflection properties of any point on the surface are those of a plane passing through that 
point, parallel to the local tangent to the surface. This method leads to a restriction on the 
rate of change of the surface gradient but no direct restriction on the magnitude of the 
surface height or surface gradient (Ogilvy 1987). 

A number of studies have considered the mechanisms governing interactions 
between underwater acoustics and rough sea surfaces. When acoustic waves interact with 
a perfectly flat sea surface, they reflect from that surface just as light would reflect from a 
mirror. That is to say, all reflected energy is directed in the specular direction. As sea 
surfaces roughen, the acoustic waves interact with the varying scales of surface 
roughness, and as a result scatter in a multitude of directions. Acoustic energy is thereby 
spread into non-specular reflections. The exact reflection parameters are dependent on the 
relationship between the length scales of surface roughness and the wavelengths of the 
impacting acoustic waves (Ogilvy 1987). 

The exact solution of acoustic waves scattering from a rough sea surface can be 
determined if the incident acoustic wave is considered to be a plane wave, while the 
surface is assumed to be infinitely periodic (Holfford 1981). 

Real ocean environments consist of a multitude of surface wavelengths and are 

not typical of these idealized environments. Therefore, a variety of scattering models 

have been utilized in the literature to approximate the interactions between acoustic 

waves and rough sea surfaces. Historically speaking, the Kirchhoff and perturbation 

methods have been the most prevalent. The Kirchhoff method is used whenever the 

acoustic wavelength is smaller than the radius of curvature of the surface. This limits the 

Kirchhoff approximation to being primarily applied to high acoustic frequencies. Because 

the Kirchhoff formulation treats all points of the scattered field on the rough surface as 

being locally flat, the approximation provides a series of reflections from different points 

along the surface. The perturbation method assumes the surface is only slightly rough, 
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and treats the rough surface as a group of perturbations from a smooth surface. As such, it 
is limited by the allowable magnitudes of surface displacement and surface slope. 
Whereas the Kirchhoff method is used primarily for high acoustic frequencies, the 
perturbation method is used mostly for low acoustic frequencies, whose wavelengths are 
generally large in comparison with the surface height (Ogilvy 1987). 

More recently, the small-slope approximation has gained in prominence. The 
small slope approximation was analyzed by Thorsos and Broschat (1995) as a series in 
generalized surface slope for the pressure-release (Dirichlet) boundary condition. This 
method was later shown, in its lowest order, to reduce to the perturbation or Kirchhoff 
approximations within each of their respective regimes (Thorsos and Broschat 1997). 

To predict the impulse responses of high frequency acoustics, rough surface 
scattering is most commonly treated through the application of the previously mentioned 
scattering models to acoustic ray codes. The main benefit of these ray codes is that they 
provide a faster calculation scheme than full-field models, such as PE solutions. A lot of 
work has been done in order to examine the effectiveness of a few surface scattering 
mechanisms applied to ray approximation models in predicting acoustic returns and 
found that while the small slope approximation gave the best results, several of the 
methods yielded inaccurate results when surface roughness increased (Williams, Thorsos, 
and Elam 2004). 

PE models, as an alternative to ray-based models, have long been used by the 
underwater acoustics community to study sound propagation. By maintaining full-field 
characteristics, PE models are able to track both the real and imaginary amplitude and 
phase values of the wave field over the entire water column, giving a total field at all 
depths and ranges. PE models are also able to track fluctuations in the acoustic channel 
over short time steps, which enable the calculation of coherence between successive 
model runs. Implementing rough surfaces into PE models follows a different approach 
than the ray codes. Rather than calculate the scattering statistics of the rough surface and 
then use that information to alter the acoustic wave, PE models instead adapt the range¬ 
marching scheme to incorporate rough surfaces into the acoustic field (Senne 2012). 
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The use of PE models for acoustic frequencies that are typical of communications 
signals has historically been sparse due to the computational effort required. However, 
recent advances in computer systems have allowed for faster calculations that enable PE 
methods to be a viable option even at moderate to high acoustic frequencies (Senne 
2012 ). 

B. DERIVATION OF PE EQUATIONS FOR THE 3-D SCATTERING 
MODEL 

Following the development of Smith (2012), we begin with a 3-D version of the 
standard PE approximation in Cartesian coordinates, 

it = ■ z>n(.*,y) (i) 


where T* is the PE “field function,” and T and U are the “kinetic” and “potential” 

op op 1 

energy operators of the PE Hamiltonian. Using the Standard Parabolic Equation (SPE) 
approximation (e.g., Smith 2001), these operators are defined in three dimensions as 


T =T 

1 op 1 3 DSPE 


1 d 2 


1 d 2 


2k 2 dz 2 2k; dy 2 


and 


U op U 3DSPE 


= ~\ (n2 


■ 1 ). 


( 2 ) 


The rough surface boundary condition is introduced by defining the location of 
the pressure release rough surface at depth 

z = rj(x,y). (3) 

The pressure release condition then requires 

X ¥(z = ?](x,y),x,y) = 0. (4) 

Introducing an image field above the surface produces odd symmetry in the field 
about the pressure release boundary which must now satisfy 

Tfi-z + 2/ j(x, y ), x, y) = -'F(z, x, y). (5) 
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Thus, the image ocean must satisfy the SPE of the form 


fi*P drj S'P 

“7 ^ ^ _ _ “z — ~rk 0 (T iDSPE + U iDSPE )'P , z<r/(x,y). 

fix ox dz 

Tappert and Nghiem-Phu (1985) introduced a field transformation 
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In the real ocean, this transformation produces 

7 )¥ , z > r/(x,y). 
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— = ~ik (T 

dx ° 3 
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( 8 ) 


In the image ocean, we must first evaluate the left-hand side of that image ocean 
form mentioned previously in Eq. (6) 
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The right-hand side of Eq. (6) includes the operator 

1 fi 2 1 fi 2 
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After combining the right-hand side with the left-hand side, we get 


dx 
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+ — 


d 2 r/ 

dxdy 


(z-rj)- 


dr/ dij 
dx dy 


dW 

dy 



dxdy 2 dxdy dy 


dri d 2 rj 
dx dy 2 


-ilk 


o 


r d 2 r/ 
y dxdy 


(z-rj) 


dr) drj^\ 
dx dy y 




( 12 ) 


Consequently, we can combine these expressions together into a single 
formulation, 

aq/ ... 

-= -ik (T +U y¥ (13) 

Qx° opop V 

with 


and 


T = 

op 


3 DSPE 


3DSPE 


il 

+- 

k 


d 2 ?] 

dxdy 


(z-rj) 


drj drj d 
dx dy dy 


, z > rj(x, y) 
, z < ? 7 (x,y) 


(14) 


u op =\ 


U 


3 DSPE 


, z > r/(x, y) 


U 


n d 2 r/, \ 

3 DSPE 2 2 (z rjj 


(15) 


„ . d } ?] d~r/ dt] drj d^i] 

+ — I- {z-rj)-2 - 


k d~i] 


dxdy 2 


dxdy dy dx dy 


dti dr] 


(z-rj)-- 

dxdy dx dy 


z < rj(x, y) 
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The previously mentioned equations represent an exact formulation of the 3-D 
rough surface scattering utilizing the SPE approximation, as originally derived by Smith 
(2012). The form of the 3-D kinetic energy operator introduces some computational 
complexities. Specifically, the field above the interface undergoes a separate cross-range 
operator. And while it is possible to apply this operator, it is not possible to apply it at a 
single depth mesh point since the rough surface varies with cross-range (i.e., q(x,y m ) ^ 
0(A,yn) )• 

In order to overcome this difficulty, an approach has been suggested by Smith 
(2013). First, at each range step (x,), the 3-D field is transfonned from cross-range spatial 
variable (y) to cross-range wavenumber variable (k y ) over all depths. Concurrently, 
transformation from depth variable (z) to vertical wavenumber variable (k z ) should also 
be performed. Second, a copy of the field in (k y ,k z ) space should be implemented. Third, 
to one copy of the field, the T op for z > q should be applied, while to the other copy the 

f for z < i] should be applied. Fourth, each copy is transfonned back to the spatial (y, z) 

domain. Fifth, based on the local surface height at range (x,), and for each cross-range 
(y m ), the appropriate form of the field should be recombined. Finally, the appropriate 
form of U should then be applied to complete the range step to (x,+i). 

The preceding 3-D expressions were developed using the SPE approximation in 
Cartesian coordinates (x, y, z). The development of the equations in cylindrical 
coordinates (r, 0, (j)) closely follows the previous approach, beginning with 

,z>7,(r,f), (16) 

where 


T =T 

1 op 1 3 DSPE 


1 d 


1 


2k; dz £ 2 k;r 2 d</> 


, U =u 

’ op 


3 DSPE 2 ^ 


- 1 ). 


(17) 
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This time the rough surface boundary condition is introduced by defining the 
location of the pressure release rough surface at depth 

z = r/(r,<f>). (18) 

Then the pressure release boundary condition is defined by 

x ¥(z = ri(r,(f>),r,<p) = 0. (19) 


Introducing an image field for z<rj produces the symmetry defined by 
'T (-z + 2 r/(r, </>), r, </>) = -'¥(z, r, <f>). 

Thus, the image ocean must satisfy a PE of the form 
S'P 3/7 <3 V P 

—+ 2 JL =- iko (T op +U op )T , z</ 7 (r,*). 

or or oz 

As before, we introduce a field transformation of the form 


( 20 ) 


( 21 ) 


x i J (z,r,<f>) = 


'h( Z,r,<f>) 


dr x ¥(z,r ,</>) 


, z > rj(r,<f>) 
, z < r/(r,<f>). 


( 22 ) 
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These equations follow exactly the same development as before, but now we end 
up with the combined formulation 


with 


~ ~ 

-= -ik (T +U ) X F 

Q r op op' 


and 


T o P =< 


1 IDS PE 


T + 

± 3DSPE 


a 


k r 2 


d~n . . dn dn 

-— (z-tj) - 

drdcj) dr d(j) 


_d_ 

d(j) 


,z> 

, z < 7(r,^) 


(23) 


(24) 


U op =\ 


u 


3 DSPE 


z > v(r,<f>) 


U 


3 DSPE 


+ - 


kj 




8 rj d 7j drj drj d rj 

-(z-rj)-2 - --i2k 


(25) 


drdfr 


drdcj) dcj) dr dcj)" 


^ d~rj drj dr] '' 

(z-t]) - 


drdcj) 


dr dcj) 


,z< ij(r,cfj) 


The implementation of these equations follows the same outline as provided in the 
Cartesian coordinates approach. 

However, it has been shown (Smith and Colosi 2008) that the SPE approximation 
performs poorly in 3-D coordinates since cross-ranges separated from the center 
introduce lateral angles that are not well treated. They suggest that a better approach in 3- 
D is to utilize the 2-D wide-angle parabolic equation (WAPE) approximation in order to 
produce the proper 3-D structure of the field. It is shown that several problems in sound 
propagation cannot be studied by using a small-angle parabolic equation, and need a 
WAPE for their solution (Ostashev and Blanc-Benon 1997). 
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Smith (2012) discussed the approach to these scattering corrections using a wide- 
angle PE. He noted that the same scattering corrections result if third order derivatives of 
the field and surface slope are ignored. Therefore, we shall utilize the 3-D WAPE 
propagator functions with the scattering correction terms derived previously. These 
operators are defined in cylindrical coordinates by 

’ u " r = u,DWArc = ~ (n ” 1) ' <26) 

Meanwhile, the WAPE introduces new terms which depend upon vertical 
derivatives of the field of second order and higher and/or products with terms that are 
second order and higher factors of the surface slope. Even though the WAPE 
implementation in the water column allows for better accuracy of high angle propagation, 
the existing application of the rough surface scattering approach does not appear to 
support the same level of high angle scattering and assumes small surface slopes. 

An approximation is done to the U in the image ocean by dropping the third 

order and higher terms in the derivative of the surface displacement that neglects the 
effects of very high slopes in the scattering and the third order and higher azimuthal 
derivatives that invoke a small-angle approximation. However, second order derivative 
tenns are maintained. This is consistent with second order derivatives of the propagating 
field, although such terms are not difficult to include. 
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III. NUMERICAL IMPLEMENTATION 


A. HYBRID IMPLEMENTATION USING A TRI-DIAGONAL SOLUTION 
FOR THE SCATTERING TERM 

The forms of the derived equations using cylindrical or Cartesian coordinates are 
identical with the obvious substitution 


d Id 

-^- . 

dy r d(p 


(27) 


For conciseness, we shall only consider the Cartesian fonn in what follows, 
although the model implementation is done using cylindrical coordinates. 

The derived equations in the previous chapter contained the additional term in the 
'/’ operator for z< //(x,yjwhich cannot be applied in the ^ wavenumber domain, as is 

traditionally done for T . Specifically, if we write 


T ={ 

op ' 


7i 


7 ; 


3DWAPE 


d 

+ R — 

dy 


, z > ri(x,y) 
, z < rj{x, y) , 


(28) 


we see that T 3DWAPE = T 3DWAPE (—-—-), which can be applied in the (k y ,k z ) wavenumber 

dy dz~ 

d d 

domain. However, R — = R(x, y,z )— cannot be treated so easily. 

dy ' dy 

d d d 

It is the combined dependence of R — on both y and —(or ^?and —in the 

dy dy d(p 

cylindrical coordinates) that does not allow a simple implementation via a Split-Step 
Fourier (SSF) algorithm. Instead, the following hybrid approach is proposed. 
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Loosely following the approach of Yevick and Thomson (1997), we will focus on 
the image field evolution equation 


a'P ~ ~ 

-= -ik H ¥ 

ax 0 0,1 


where, H is a Hamiltonian-like operator that defines the evolution of the PE field 
function in range. We may then use the approximate marching algorithm solution 


(29) 


Tfix + Ax) = e 


-ik A rH (r) 


'P(x). 


(30) 


Here, 

H . r w = L w+ u ,w = *+ K + u* • O'i) 

so we may define, 

_ -ik fJ AxR— _ 

Y(x + Ax) = e ^ e -'K^ De ~,KMU, DX¥{x) ' ( 32 ) 

Note that this form assumes that the rough surface correction is applied after the 
kinetic energy operator (T 3D ). 


The Monterey-Miami Parabolic Equation (MMPE) model employs a centered- 
step scheme implementation, which in the real ocean takes the form (Smith and Tappert 
1993, Smith 2001) 


-i—k„AxU-, n (x+Ax) , 

X ¥(x + Ax) = e 2 FFT;! 


-ik Q AxT 3D ( k x ,k y ) 


FFT. 


— i - k \xlI -) r\ 

e 2 'P(x) 


(33) 


where 'P(x) = 'P(x,y,z)and T 3D (k x ,k y ) is the scalar form of the operator T, L) in the (k y , 
k z ) wavenumber space. 
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The image ocean solution can utilize the SSF algorithm for the T 3D term, but must 


employ an alternative method for the R —term. Thus, we will implement a centered-step 

dy 

scheme defined by 


-i-k 0 AxU 1D (x+Ax) -ik a AxR(x+]-Ax)-^- 


X ¥(x + Ax) = e 


2 dy J7JJ'J 1 J (k x ,k y ) 


FFT 


-i-k o AxU, D (x) 


T(x) 


>• (34) 


In order to implement this equation, we should determine how to apply the 


operator e ,?Rr with g = k Ax and y = — by defining the intermediate step solution as 

dy 


T'U + Ax') = FFT;: FFT 


~K;SU, D (x) 


'T(.x) 


(35) 


The final solution at the end of each range step is then obtained from 


~ -i—cU-> n (x+Ax) ~ 

X ¥(x + Ax) = e 2 T / (.x + Ax ) , 


(36) 


where 


V(x + Ax') = e- iiRr V'(x + Ax'). 


(37) 
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The rough surface correction operator can be represented by the Pade 
approximant (Collins 1990) 


e - igRy « 


i ~\^ R y 


i + ^igRy 


(38) 


which is valid for small y , implying small corrections in cross-range scattering at each 
step. From this approximation, we obtain a finite difference system of equations based 
upon 


1 


1 

l + -ig R r 

X ¥(x + Ax') = 



'F(jc + Ax'). 


This may be written in a vector-matrix form as 


A^ = BW'=C , 


where 


"V 



'P 


V P' 

2 


2 


T" = 






and ^ ='P(y n ) , n=l,2,3....N. 


(39) 


(40) 


(41) 
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The expressions required for the A and B matrices’ components will be based on a 
3-point centered approximation to the derivative 


dy 


1 


y o 


2 Ay 




(42) 


where v F +] = T(y = y ± Ay). At the end points (n=l or n=N), we shall invoke the 2-point 
approximations 


f¥ = 


5^ 

dy 


y a 


A? 




dy 


y a 


— (A' -T ). 
Ay 1 1 ° } 


The former is used if n=N and the latter applies when n=l. 
The matrices may be defined in terms of a and b where 

J n n 


i + ^igRr 


T ='F + -icRW = 

O O 2 ' O 


'P +-igR— (T,-T ,) = ¥ +-^-('P,-T ,), 
0 2 2Ay 1 - l 0 4Ay 1 - 1 


from which we can introduce 


(43) 


(44) 


(45) 


1 'SK 

a - — 


" 4 Ay 


b = 


"I 'S R n 

4 Ay 


R =R(y ). 

n n 7 


(46) 
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Consequently, we can define two different matrices 


and 


" (1-2a,) 2a, 0 0 0 0 

-a 2 1 a 2 0 0 0 

0 -a 3 1 a 3 0 0 

0 0 0 0 0 
0 0 0 -<V> 1 a N-\ 

v 0 0 0 0 -2 a N (1 + 2 a N ) J 


" ( 1 - 2 ^) 2 b x 0 0 

-b 2 1 b 2 0 

0 -b 3 1 b 3 

0 0 0 

0 0 0 -b N _ x 

0 0 0 0 


0 

0 

0 

0 

1 


0 

0 

0 

0 

V i 

(1 + 2 **) , 


Now the matrix equation becomes 


A'F = C = 


C 3 


( 47 ) 


(48) 


(49) 


C 


N 


with the elements of C are computed from BW ■ Equation (50) can be solved by a tri¬ 
diagonal matrix algorithm (Press et al. 2007) according to the recursive relations 


'P =C' 

W .=a-aW. 


i +1 


i = 1,2,3. 


,N-2,N-1 , 


(50) 
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where 



2 a x 

1-2 ~a x 
a 

_ l _ 

1 + a' .a. 

I-l l 


i=l 


i=2,3,...,N-l 


C' = 


C, 


1-2 a, 


C. + C a. 

i _ / i i 

1 + a' ,a. 

l-l l 


C N + 2 ci C' , 

_ N _ N N-\ 

(\ + 2a N ) + 2a N a' N1 


i=l 


i=2,3,...,N-l 


i=N . 


( 51 ) 


(52) 


The outcome of this tri-diagonal matrix solution is then x F(x + Ax') to which we 

-i-5U ?> DSPE^ x+ A*) 

simply apply the scalar vector multiplier e 2 to obtain the solution at the end of 

the range step ^(x + Ax). 
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B. DEFINITION OF ROUGH SURFACE PERTURBATIONS 


In order to appropriately test the 3-D scattering model just developed, it is 
required that we analytically define some simple surface perturbations. In this section, we 
shall define five simple surface perturbations. These will be defined in cylindrical 
coordinates, since the model is implemented this way. 

First, we define a surface with constant slope in order to determine the way the 
different radials will be affected by that surface slope using the equation 


r/(r,(/>) = 20 + 0. l(r cos (/>), (53) 

where rj(r,(/)) is the surface height with respect to range, r , and the expressions for the 
required derivatives are then 


— = O.lcos(^), 
dr 


= 0 , 


d 2 r] 
dr 2 

— = -O.lrsin(^), 
d<f> 


d 2 r/ 

drd(j) 


= -O.lsin(^). 


(54) 


This surface is depicted in Fig. 1. 

Second, we define the surface roughness for concentric sinusoidal waves 

rj{r , (f)) = <J cos (kr), (55) 


where ij(r,(/)) is now the surface wave height with respect to range, r (independent of 
bearing (j) in this case), (7 is the surface wave height amplitude, and k is the surface 
wavenumber (2 njA s where A s is the wavelength of the surface perturbation). Notice here 

that since the surface does not depend upon the bearing angle, the propagation out along 
all radials should be affected similarly. Then 
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di] 

dr 

d 2 ?] 


= -kcrsm(kr), 

, =-k 2 crcos(kr), 


dr 

*L = o, 

d(/> 


d 2 r/ 

drd<f> 

This surface is depicted in Fig. 2. 


= 0 . 


( 56 ) 
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Figure 2. Concentric sinusoidal surface perturbations. 


Third, we define the equations introducing the surface roughness for 
unidirectional sinusoidal waves by 

rj(r,(f>) = acos(krcos</>), (57) 


where (f> is the angle of propagation along each radial. Then 


—- = -kcr cos(^) sin[/:rcos(0)], 
dr 

d 2 rj 

—Y = -a[k cos(^)] 2 cos[&r cos(^)], 
dr 

dii 

— = krtj sin(^) sin [kr cos(^)], 
d(j) 

d 2 n 

- ^ = ka sin(^) |sin[/cr cos(^)] + kr cos ((f>) cos[£r cos(^)]}. 


(58) 
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Figure 3. Planar sinusoidal surface perturbations. 


Fourth, we define a single upward sinusoidal perturbation on one side only by 
using the same surface roughness equation as the unidirectional surface waves, but this 
time we limit the angle of propagation (f) to be within the forward region of the source and 
only allow one oscillation, 


7j(r, <j)) = <J cos {kr cos cp) 


-90 < (j) < 90, 
r cos (j) < X s . 


(59) 
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Elsewhere, the surface will be constant with a value for r/(r,</>) equal to the 
maximum surface wave height cr. The derivatives in the region of the displacement are 
then 


—- = -£crcos(^)sin[/:rcos(0)], 
dr 

p\1 

—% = -<j[k cos(^)] 2 cos[/crcos(^)], 

dr" 

— = &rasm(^)sin[/:rcos(0)], 
d(j) 

- = kcr sin(^) {sin[/:r cos(^)] + kr cos {(/)) cos[£r cos(^)]}, 

drd(j) 


while over the other directions and ranges, the sea surface is defined as 


r](r,(f>) = a. 


dr/ 

dr 

d 2 r/ 

dr 2 

dr/ 

d(j) 


= 0 , 
= 0 , 
= 0 , 



drd(j) 


This surface is depicted in Fig. 4. 


(60) 


(61) 


Finally, we flip the single upward sinusoidal perturbation wave and shift it until it 
reaches the mean sea level, where r/{r,(f>) equals zero everywhere except in the forward 
region of the source, where the sea surface equation is defined by 


-90 <(/)< 90, 

t/(r, < />) = cr[l-cos(krcos(t))] , (62) 

r cos tp<A s . 
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Figure 4. Upward sinusoidal surface perturbations. 


The derivatives in the region of the displacement are then 


3/7 

—- = kcr cos (0) sin [kr cos(^)], 
dr 

= a[k cos(^)] 2 cos [kr cos(^)], 
dr 

dii 

—— = -kr<j sin(^) sin[/crcos(^)], 
d(j) 

d 2 n 

^ ^ = ~k<J sin(^) |sin[/:r cos(^)] + kr cos(^) cos[£r cos(^)]}. 


while all over the other directions, the sea surface is defined as 


(63) 
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r/(r,0) = O, 


dr/ 

dr 


= 0 , 


d 2 r/ 

d7 


= o, 


dr/ 

d(j) 


= 0, 


d 2 r/ 

drdrj) 

This surface is depicted in Fig. 5. 


(64) 
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IV. MODEL RESULTS 


A. CONSTANT SLOPE SURFACE VERSUS CONSTANT SLOPE BOTTOM 

Although a sloping surface is a relatively benign definition of a perturbed surface, 
it does allow a comparison of the rough surface algorithm defined here with a formally 
equivalent environment defined by a sloping bottom. The solution to a constant slope 
surface is identical to a constant slope bottom through a simple environmental rotation, 
assuming ideal point sources in both cases. The approach to the 3-D problem with a non¬ 
constant bottom is no different from the approach to the 2-D problem, where the MMPE 
model simply adjusts the bottom depth at each range step and radial. 

1. Constant Slope Surface and Flat Bottom 

Solutions to the problem of a sloping surface over a flat bottom are depicted in the 
following figures. The source is located at a depth of 50 m and the bottom depth is 
highlighted with a white line at 120 m depth. Figures 6 and 7 display the solutions along 
opposing radials, with the former along the direction of descending surface depth and the 
latter along the direction of ascending surface depth. In both cases, the bottom interface is 
observed to be flat. 

In Figure 8, the transmission loss at a single depth, corresponding to the source 
depth, is extracted from all radials. Note that in the case of a perfectly flat surface and 
bottom, this would produce concentric rings of interference patterns in the field. In this 
case, however, the directionality of the propagation along different radials is apparent. 
Transforming this data to a polar coordinate system produces the results presented in 
Figure 9. The effect observed is consistent with general expectations of a sloping surface, 
where the propagation in the direction of the descending surface produces interference 
patterns compressed in range and, alternatively, the propagation in the direction of the 
ascending surface produces interference pahems elongated in range. 
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Depth (m) Depth (m) 


Transmission Loss (dB re 1m) 



Range (km) 

Figure 6. 0-degree radial for a source at depth of 50 m. 


Transmission Loss (dB re 1m) 



Range (km) 

Figure 7. 180-degree radial for a source at depth of 50 m. 
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Transmission Loss (dB re 1m) 
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Figure 8. Transmission loss pattern at source depth. 


Transmission Loss (dB re 1m) 



Range (Km) 

Figure 9. Constant slope surface effect on the waves emitted by a 50 m deep source. 
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2. Constant Slope Bottom and a Flat Surface 

The corresponding solution is now computed by interchanging the surface and 
bottom characteristics. The environment now has a smooth flat surface with a sloped 
bottom which should produce the equivalent field solution, albeit rotated, as in the 
constant slope surface case. In this case, the surface depth is set to a constant value of 
zero and the source depth is shifted upwards to 30 m depth (consistent with source depth 
relative to the surface in the previous case), while the average depth of the medium will 
be 100 m. The 2-D bottom bathymetry for this environment is defined as 

bdint = bdint - 0. l(rcos^) (65) 

where “bdint” is the interpolated bottom depth at each range step and (j> is the bearing 
angle of propagation. 

Figures 10 and 11 present the results along opposing radials. The results in these 
figures should compare directly, through a simple rotational transformation, to the results 
presented earlier in Figures 6 and 7. 
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Figure 10. O-degree radial for a source at depth of 30 m. 
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Figure 11. 180-degree radial for a source at depth of 30 m. 

3. Constant Slope Surface Model versus Constant Slope Bottom Model 

In order to compare results from the two different approaches, a common depth 
for extracting transmission loss must be defined. The most convenient depth to choose in 
both cases is the depth of the bottom bathymetry. Other than some geometrical rescaling 
(due to rotation) of the range from the source, the transmission loss traces along the 
bottom depth should agree between these two methods. 

Figure 12 shows the initial results directly extracted from each solution at the 
bottom interface in polar coordinates. Due to the inclination presented by the sloping 
bottom model, we must adjust the range defined in the calculation in order to compare 
results. Figures 13 and 14 utilize this range-scaling correction to compare the 
transmission loss traces along opposing radials. 

As can be noted in Figures 13 and 14, the comparison between solutions 
computed with a 2-D sloping bottom using the standard 3-D MMPE and this new 
algorithm with a 2-D sloping surface show good agreement. 
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Figure 12. 


Model solutions for the water-bottom interface: (left) the constant slope surface; 

(right) the constant slope bottom. 



Range(Km) 


Figure 13. 


Transmission loss at water-bottom interface along the 0 radial. 
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Range(Km) 

Figure 14. Transmission loss at water-bottom interface along the 180 radial. 

B. CONCENTRIC SINUSOIDAL SURFACE PERTURBATIONS 

As an alternative approach to confirm proper model behavior, we now examine 
the results derived from a surface perturbation defined as concentric sinusoidal waves. 
Theoretically, all the radials should be affected in the same manner due to the defined, 
radially-independent environment. 

Figures 15 and 16 present the transmission loss field results from opposing 
radials. As expected, the results appear identical in their structure along each radial. 
Figure 17 displays the transmission loss results extracted from a constant depth of 30 m, 
the same as the source depth. The results transformed to polar coordinates are presented 
in Figure 18. 

The results show, as expected, that the acoustic field will propagate identically 
along all bearing angles, producing the circular patterns shown in Figure 18. 
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Figure 15. O-degree radial for a source at depth of 30 m. 
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Figure 16. 180-degree radial for a source at depth of 30 m. 
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Figure 17. Transmission loss pattern at source depth. 
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Figure 18. Concentric sinusoidal surface perturbations effect on a transmission loss field 

extracted at source depth. 
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c. 


PLANAR SINUSOIDAL SURFACE PERTURBATIONS 


We now examine the results computing using the planar sinusoidal surface 
perturbations. While the forward and backward (0 and 180 degrees) radials should 
produce the same results due to symmetry, there should be significant differences 
between those and other radials, where the surface perturbation along the path is 
different. 
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Figure 19. O-degree radial for a source at depth of 30 m. 
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Figure 20. 180-degree radial for a source at depth of 30 m. 
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Figures 19 and 20 display the symmetric radials along 0 and 180 degrees. As 
expected, the transmission loss patterns appear identical between the two bearings. 
However, since the propagation will respond differently along each radial, an 
examination of the field at a single depth produces the results displayed in Figures 21 
and 22. 

Note that the dominant scattering features present (most notably in Figure 22) 
occur upon the first interaction with the downward perturbation of the surface roughness. 
The impact on the field is clearly parallel to the plane of the surface perturbation, and for 
ranges shorter than +/- 0.1 km on the horizontal cross-range axis, energy appears to be 
channeled between the downward portions of the surface roughness on either side. 

Beyond +/- 0.1 km on the horizontal cross-range axis, the impact of the surface 
roughness is less noticeable. This is due to the dominant scattering occurring within the 
first +/- 0.1 km. Beyond this range, much of the surface interaction has been shadowed by 
the near field scattering. 
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Figure 21. Transmission loss pattern at source depth. 
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Figure 22. Planar sinusoidal surface perturbations effect on a transmission loss field 

extracted at source depth. 


D. UPWARD SINUSOIDAL SURFACE PERTURBATIONS 

In order to isolate some of the scattering phenomena occurring from the previous 
type of planar sinusoidal surface perturbation, and to break the symmetry of the problem 
to ensure proper model behavior, the next surface to test is simply a single sinusoidal 
surface displacement. We start initially by creating a single upward sinusoidal surface 
perturbation to the field, identical in form to the previous surface on one side only and 
out to 0.1 km. We select the 0 radial to represent the source front-side of the acoustic 
propagation and the 180 radial to show the response on the back-side of the source. 

The results are presented in Figures 23 through 26. As expected, the back-side of 
the propagation shows little impact from the surface perturbation on the front-side. There 
may be some minor effect along the boundary separating front-side and back-side, but 
this is to be expected due to horizontal scattering. In addition, the scattering in the front¬ 
side direction is very similar to the previous multiple sinusoidal planar perturbation, 
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confirming that most of the scattering effects are due to the initial interactions with the 
downward sloping portion of the surface. 



Figure 23. 0-degree radial for a source at depth of 30 m. 


Transmission Loss (dB re 1m) 

10 
20 
30 
40 
50 
60 
70 
80 
90 
100 

0.05 0.1 0.15 0.2 0.25 0.3 

Range (km) 

Figure 24. 180-degree radial for a source at depth of 30 m. 
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Figure 25. Transmission loss pattern at source depth of 30 m. 
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Figure 26. Upward sinusoidal surface perturbations effect on a transmission loss field 

extracted at source depth. 
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E. DOWNWARD SINUSOIDAL SURFACE PERTURBATIONS 

As an alternative to the single sinusoidal planar perturbation presented previously, 
we also created a similar surface that perturbs the surface downward initially before 
returning upward to a constant depth. One of the primary differences of this perturbation 
was that it places the strongly scattering downward slope of the surface perturbation 
immediately adjacent to the source range. This would, therefore, be expected to introduce 
stronger scattering at higher angles, particularly in the near-field. As before, we choose 
the 0 radial to represent the front-side path of the acoustic propagation and the 180 radial 
to represent the back-side field of the source. 

It is worth noting that during early testing of the model’s stability, this particular 
surface perturbation initially generated non-physical results. Closer inspection revealed 
that the earlier version of the model was utilizing an SPE approximation to compute the 
azimuthal wavenumber propagator function, resulting in large scattering angles being 
misrepresented in the solution. As discussed in the theory section of this thesis, it is quite 
important to employ wide-angle approximations, even in cross-range/bearing-angle, in 
order to properly compute the azimuthal field structure. After upgrading the code to 
account for this, the results from this perturbation model became physically realizable. 
(All results reported in this thesis, therefore, utilized the WAPE approximation for all 
operators.) 

Figures 27 through 30 display the computed transmission loss field results for this 
surface perturbation. Again, the back-side of the field shows little impact from the 
perturbation, except very near the edge separating the front-side from the back-side, 
which is to be expected. In addition, the front-side scattering is much stronger in this case 
than in the previous one. There is also clearly diffraction of the field filling in the 
classical shadow zone on the back-side of the downward surface perturbation. 

It is natural to expect that the scattering from the front-side would eventually be 
transferred across the boundary to the back-side of the field. Because of the limitation on 
the source angle in the vertical, and the relatively short range of the calculations, this was 
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not clearly observed in any of these calculations. It is worth exploring this issue further in 
future analysis. 
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Figure 27. 0-degree radial for a source at depth of 30 m. 
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Figure 28. 180-degree radial for a source at depth of 30 m. 
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Figure 29. Transmission loss pattern at source depth. 


Transmission Loss (dB re 1m) 



Range (Km) 


Figure 30. Downward sinusoidal surface perturbations effect on a transmission loss field 

extracted at source depth. 
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F. 


COMPARISON OF THE EFFECT OF 3-D VERSUS N x 2-D 


One of the goals of this project is to develop an accurate means of computing the 
3-D scattering from a rough 2-D surface, and then comparing solutions between this 3-D 
model and standard Nx2-D models. Such Nx2-D models are significantly more efficient 
and require much less in the way of computing power. And so, understanding when a true 
3-D solution is required will help avoid unnecessary computational complexity. In this 
section, we will compare results computed between an Nx2-D version of the MMPE 
model that utilizes a consistent rough surface scattering algorithm with the 3-D version 
described in this thesis. 

1. Constant Slope Surface Perturbations 

Utilizing the same environment defined previously for the constant sloping 
surface, the results from the 3-D model and the corresponding Nx2-D version were 
computed out to a maximum range of 0.6 km. The goal of carrying out these calculations 
to a slightly longer range than before is simply to allow for more azimuthal coupling to 
occur which would differentiate the impact of the 3-D calculation. Figures 31 through 34 
present the model results from this environment along the different radials and at a 
constant depth. 

The results from the two models, particularly as displayed in Figure 34 in polar 
coordinates, show little difference between the models. This is perhaps not too surprising, 
since azimuthal coupling is known to have a larger impact at longer ranges. More 
significant differences would likely be observed if the calculations were carried out at a 
lower frequency and the maximum range was increased significantly. However, due to 
the nature of the sloping surface, longer ranges would cause the surface to cross the 
bottom interface, generating nonsensical environments. 

In an effort to highlight the differences computed from the 3-D and Nx2-D 
models, Figure 35 displays the difference in transmission loss values at a constant depth. 
As noted previously, the differences are relatively minor, and only begin to become 
prominent at longer ranges. Furthermore, the strongest differences occur along the radials 
at 90 and 270 degrees, corresponding to bearings of constant surface depth. These 
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bearings also correspond to the maximum cross-range slopes of the surface, where we 
would expect to see the strongest amount of azimuthal coupling. 
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Figure 31. Model solutions for the constant slope surface perturbations along the 0 radial: 

(left) using the basic 2D interface perturbations model; (right) using the 3D 

model. 
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Figure 32. Model solutions for the constant slope surface perturbations along the 180 radial: 

(left) using the basic 2D interface perturbations model; (right) using the 3D 

model. 
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Figure 33. The wave pattern results for a 50 m deep source: (left) using the basic 2D 
interface perturbations model; (right) using the 3D model. 
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Figure 34. Model solutions for the constant slope surface perturbations effect on a 50 m deep 
source waves: (left) using the basic 2D interface perturbations model; (right) 

using the 3D model. 
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Figure 35. Difference in transmission loss values at a constant depth. 

2. Downward Sinusoidal Surface Perturbations 

The second environment we will consider for comparison between the fully 3-D 
results versus the Nx2-D approach corresponds to the single, downward sinusoidal 
surface perturbation. This perturbation produced the strongest scattering in the near-field, 
and would be expected to produce strong differences between 3-D and Nx-2D in the 
near-field, as well. 

Figures 36 through 39 present the results from the two approaches for this 
environment. It is worth noting that the single radial results presented in Figure 36 show 
some marked differences in the field computed in the bottom when the full 3-D model is 
employed. This is due to high vertical angle scattering produced by the computational 
algorithm. The transmission loss levels of these differences are very low, however, nearly 
120 dB below the primary energy observed in the water column, and are not of any 
physical significance. 
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Figure 36. Model solutions for the downward sinusoidal surface perturbations along the 0 
radial: (left) using the basic 2D interface perturbations model; (right) using the 3D 

model. 
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Figure 37. Model solutions for the downward sinusoidal surface perturbations along the 180 
radial: (left) using the basic 2D interface perturbations model; (right) using the 3D 

model. 


As before, the polar plots of the transmission loss field at a constant depth of 30 m 
displayed in Figure 39 shows little apparent difference in the solutions. In an effort to 
highlight the differences, Figure 40 shows the difference in transmission loss at a 
constant depth in polar coordinates. As before, the differences are observed to expand in 
range, and are strongest in regions where the cross-range slope of the surface perturbation 
is largest. 
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Figure 38. The wave pattern results for a 30 m deep source: (left) using the basic 2D 
interface perturbations model; (right) using the 3D model. 
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Figure 39. Model solutions for the downward sinusoidal surface perturbations on a 30 m 
depth source waves: (left) using the basic 2D interface perturbations model; 

(right) using the 3D model. 
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Figure 40. Difference in transmission loss values at a constant depth. 
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V. CONCLUSIONS AND FUTURE WORK 


Rough surface scattering is an important mechanism for introducing signal 
variability into underwater acoustic systems such as acoustic communications. The 
significance of three-dimensional scattering effects from two-dimensional rough surfaces 
has been investigated by other researchers. However, whether the full 3-D calculations 
constitute a significant improvement to the simpler, Nx2-D calculations remains unclear. 
The goal of this thesis was then to upgrade an existing 2-D parabolic equation model to 
incorporate full 3-D surface scattering. This will allow future studies that can more 
directly compute the specific influence of 3-D, out-of-plane scattering from 2-D rough 
surfaces. 

There were two primary advances made in this work towards the stated goal. 
First, a hybrid implementation for the parabolic equation model was formulated and 
incorporated into a version of the Monterey-Miami Parabolic Equation model. This 
hybrid method consisted of the standard split-step Fourier algorithm traditionally 
employed in MMPE, and also a finite difference method utilizing a Pade approximation 
for the cross-range operator associated with the cross-range scattering of the rough 
surface. The solution to the finite difference system of equations was accomplished 
through a standard tri-diagonal solver. 

The second advance related to the wavenumber propagator functions employed by 
the 3-D model in cylindrical coordinates. Previous work assumed the azimuthal coupling 
was small enough that a standard parabolic equation approximation was sufficient. 
However, with some of the larger surface slopes investigated in this thesis, it was 
determined that a wide angle parabolic equation approximation was necessary in order to 
properly account for the cross-range scattering. Once this upgraded propagator was 
incorporated, the model results became physically realistic for all of the rough surfaces 
investigated. 

Five different “rough” surfaces were defined for model testing. These were all 
simple, non-flat boundaries that increased in complexity. The first example studied was 
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simply a smooth, sloped surface. The resulting field appeared consistent with 
expectations. In order to validate these results, the model was recast to compute the 
solution for a flat surface with a sloping bottom. Through an environmental rotation, 
these results should be identical. A comparison of the model results confirmed that the 
sloping surface results produced valid solutions. 

The second surface was defined by concentric, sinusoidal surface displacements. 
Being independent of bearing angle, this surface was designed to test the model 
calculations and ensure that no erroneous out-of-plane scattering was occurring when 
none was expected. The model results were consistent with this expectation. 

The third, fourth, and fifth surface fluctuations were defined in terms of 
horizontally planar sinusoidal surface perturbations of varying types. While there was no 
formal solution to compare with the model results, the general behavior was consistent 
with expectations, where scattering was observed parallel to the surface wavefront. 

Finally, model results for two of the surface perturbations were computed using 
both the 3-D version developed here, as well as the Nx2-D version previously developed. 
A comparison of the results showed the dominant differences along the directions with 
the greatest cross-range surface slopes, as expected. 

Future work could entail longer range solutions in order to search for more 
obvious cross-range coupling. Using sources defined by individual waveguide modes 
could also be employed to more clearly delineate horizontal refraction. In addition, the 
model results examined in this thesis were computed using only a single frequency. 
Future work could examine broadband signals, and thereby compute pulse propagation 
travel times, which could also provide specific evidence of horizontal refraction and out- 
of-plane scattering. 
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